rm(list = ls())

#设置工作路径
setwd("Replication_BRI")

library(cshapes)
library(splm)
library(dplyr)
library(purrr)
library(ggplot2)
library(ggthemes)
library(tidyr)
library(plm)
#加载自制函数用于可视化
source("0_setup.R")
#加载数据
load("data.RData")

## 模型1-4：图2-图3

m00 <- glm(bri_join_dum ~ initiative_cnty + neighjoin00_bri_lag1+ 
             years_diplomatic_lag1 +  IdealPointDistance_chn_lag1 + 
             arms_logexportfromCHN_lag1+  arms_logimporttoCHN_lag1 +
             v2x_polyarchy_lag1+  ctygdppc_log_lag1 + dep_import_lag1 +
             dep_export_lag1 + logviolences_lag1 + factor(year),
           family = binomial(link = "logit"), data = data)

m50 <- glm(bri_join_dum ~ initiative_cnty + neighjoin50_bri_lag1+ 
             years_diplomatic_lag1 +  IdealPointDistance_chn_lag1 + 
             arms_logexportfromCHN_lag1+ arms_logimporttoCHN_lag1 + 
             v2x_polyarchy_lag1+  ctygdppc_log_lag1 + dep_import_lag1 + 
             dep_export_lag1 + logviolences_lag1 + factor(year),
           family = binomial(link = "logit"), data = data)

m200 <- glm(bri_join_dum ~ initiative_cnty + neighjoin200_bri_lag1+ 
              years_diplomatic_lag1 +  IdealPointDistance_chn_lag1 + 
              arms_logexportfromCHN_lag1+ arms_logimporttoCHN_lag1 +
              v2x_polyarchy_lag1+  ctygdppc_log_lag1 + dep_import_lag1 +
              dep_export_lag1 + logviolences_lag1 + factor(year),
            family = binomial(link = "logit"), data = data)

m400 <- glm(bri_join_dum ~ initiative_cnty + neighjoin400_bri_lag1+ 
              years_diplomatic_lag1 +  IdealPointDistance_chn_lag1 + 
              arms_logexportfromCHN_lag1+ arms_logimporttoCHN_lag1 +
              v2x_polyarchy_lag1+  ctygdppc_log_lag1 + dep_import_lag1 +
              dep_export_lag1 + logviolences_lag1 + factor(year),
            family = binomial(link = "logit"), data = data)

library(texreg)
library(showtext)

# 图2（b)
logit_f_coef <- scaled_coef_cluster(ModelResults = list(m00, m50, m200, m400), 
                                    data = data, 
                                    clusterid = "iso", subvar = TRUE, 
                                    vars = c("neighjoin00_bri_lag1", 
                                             "neighjoin50_bri_lag1",
                                               "neighjoin200_bri_lag1",
                                              "neighjoin400_bri_lag1"))

logit_f_coef + scale_x_discrete(labels = parse(text = c("邻国参与数量[0~km]",
                                           "邻国参与数量[50~km]",
                                           "邻国参与数量[200~km]",
                                           "邻国参与数量[400~km]")))
ggsave("Fig2_b.png", width = 6, height = 4.2)

# 安装开发版的包
library(postregplots)

## 图-2 a
predict_probs(ModelResults =  m400, n.sim = 1000, 
              varname = "neighjoin400_bri_lag1", 
              data =data, val1 = 0,
              val2 = 18, intervals = 1, clusterid = "iso") + 
  labs(title = "",
       x = "邻国参与数量",
       y = "预测概率")+ 
  theme(plot.title = element_text(hjust = 0.5, size = 16),
        strip.text = element_text(size = 15),
        text = element_text(size=16, family = 'KaiTi',face = "bold")) 
ggsave("Fig2_a.jpeg", width = 4, height = 3)

### 图3-a

coef_comparison2(ModelResults = list(m00, m50, m200, m400), 
                 n.sim = 1000, varname = c("neighjoin00_bri_lag1", 
                                            "neighjoin50_bri_lag1",
                                             "neighjoin200_bri_lag1",
                                             "neighjoin400_bri_lag1"),
                data = data,
                val1 = 0.025, val2 = 0.975, clusterid = "iso") +
            theme(plot.title = element_text(hjust = 0.5, size = 16),
                  strip.text = element_text(size = 14),
                  text = element_text(size=14, family = 'KaiTi',face = "bold")) +
            scale_y_discrete(labels = parse(text =c( 
              "邻国参与数量[0~km]",
              "邻国参与数量[50~km]",
              "邻国参与数量[200~km]",
              "邻国参与数量[400~km]"))) + 
            labs(title = "", x = "平均边际效应") 
ggsave("Fig3_a.jpeg", width = 4, height = 4)

## 图3-b

coef_comparison(ModelResults = list(m00, m50, m200, m400), n.sim = 1000,
                varname = "initiative_cnty",
                data = data,
                val1 = 0, val2 = 1, clusterid = "iso") +
        theme(plot.title = element_text(hjust = 0.5, size = 16),
              strip.text = element_text(size = 14),
              text = element_text(size=14, family = 'KaiTi',face = "bold")) +
        scale_y_discrete(labels = parse(text =c( 
          "模型~1[0~km]",
          "模型~2[50~km]",
          "模型~3[200~km]",
          "模型~4[400~km]"))) + 
        labs(title = "", x = "平均边际效应") 
ggsave("Fig3_b.jpeg", width = 4, height = 4)


##### 附录图-1
library(ggcorrplot)
library(extrafont)
library(showtext)
library(ggcorrplot)
font_add('KaiTi', 'KaiTi.ttf')
df <- data %>% ungroup() %>% 
  dplyr::select(bri_join_dum, initiative_cnty, neighjoin00_bri_lag1, years_diplomatic_lag1 ,  
             IdealPointDistance_chn_lag1 , arms_logexportfromCHN_lag1, 
             arms_logimporttoCHN_lag1 ,v2x_polyarchy_lag1,  ctygdppc_log_lag1 ,
             dep_import_lag1 , dep_export_lag1 , logviolences_lag1)

cor_df <- na.omit(df)
corr <- round(cor(cor_df ), 1)

ggcorrplot(corr, hc.order = F, type = "lower", 
                colors = c("blue", "white", "red"),digits = 2,
            lab = TRUE, sig.level = 0.05, ggtheme = ggplot2::theme_minimal) + 
          ggplot2::theme(legend.title=element_blank(),
                axis.text = element_text(size=10),
                text = element_text(size=10,family ='KaiTi'))+
            ggplot2::scale_x_discrete(labels = parse(text = c("邻国加入数量[list(t-1)]",
                                            "建交历史长度[list(t-1)]",
                                            "与中国外交立场差距[list(log,t-1)]",
                                            "向中国进口军售额[list(log,t-1)]",
                                            "向中国出口军售额[list(log,t-1)]" ,
                                            "政体类型[list(t-1)]",
                                            "人均GDP[list(log, t-1)]",
                                            "对中国出口贸易依赖度[list(t-1)]",
                                            "对中国进口贸易依赖度[list(t-1)]",
                                            "国内暴力数量[list(log,t-1)]")))  +
        ggplot2::scale_y_discrete(labels = parse(text = c( "初始沿线国家[list(t-1)]",
                                                     "邻国加入数量[list(t-1)]",
                                                     "建交历史长度[list(t-1)]",
                                                     "与中国外交立场差距[list(log,t-1)]",
                                                     "向中国进口军售额[list(log,t-1)]",
                                                     "向中国出口军售额[list(log, t-1)]" ,
                                                     "政体类型[list(t-1)]",
                                                     "人均GDP[list(log,t-1)]",
                                                     "对中国出口贸易依赖度[list(t-1)]",
                                                     "对中国进口贸易依赖度[list(t-1)]"))) 
                g 
ggsave("corr.png", width = 10, height = 10)

# 附录表-1
library(stargazer)
df <- as.data.frame(df)
stargazer(df, header=FALSE, out = "table.txt",  type = "text", title = "描述性统计结果",
          digit.separator = "")

